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We consider propagation of surface plasmon polaritons (SPPs) in linear periodic chains (LPCs) 
of prolate and oblate metallic spheroids. We show that the SPP group velocity can be efficiently 
controlled by varying the aspect ratio of the spheroids. For sufficiently small aspect ratios, a gap 
appears in the first Brillouin zone of the chain lattice in which propagating modes do not exist. 
Depending on the SPP polarization, the gap extends to certain intervals of the Bloch wave number 
q. Thus, for transverse polarization, no propagating SPPs exist with wave numbers q such that 
Qc < \l\ < ft/h, h being the chain period. For longitudinally polarized SPPs, the gap spans the 
interval \q\ < q\. Here qjr and q\ are different constants which depend on the chain parameters, 
spheroid aspect ratio and its orientation with respect to the chain axis. The dependence of the 
dispersion curves on the spheroid aspect ratio leads to a number of interesting effects. In particular, 
bandwidth of SPPs that can propagate in an LPC can be substantially increased by utilizing prolate 
or oblate spheroids. When q is close to a critical value, so that \q — q^\ <C Tf/h or \q — qj <C ir/h, 
' the decay length of the SPPs is dramatically increased. In addition, the dispersion curves acquire 

a very large positive or negative slope. This can be used to achieve superluminal group velocity for 
realistic chain parameters. We demonstrate superluminal propagation of Gaussian wave packets in 
numerical simulations. Both theory and simulations are based on Maxwell equations with account 
of retardation and, therefore, are fully relativistic. 
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I. INTRODUCTION 
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Propagation of surface plasmon polaritons (SPPs) in linear periodic chains (LPC) of metal nanoparticle has been in 
£>V the focus of considerable recent attentioriii^i&ii^i&i^. The interest is, in part, motivated by the potential application 
*^ of such chains as sub-wavelength plasmonic waveguide o 9 i 10 ' 1 1 . Since energy in a waveguide can be transported in the 
i ^~f ' form of wave packets, the dispersion relation becomes of primary importance. In the case of LPCs, the dispersion 
relation is a mathematical dependence of the SPP frequency to on its Bloch wave number, q. It is important to 
CN| ■ emphasize that the elementary excitations (electromagnetic modes) in LPCs are Bloch waves rather than sinusoidal 
waves which propagate in continuous media. Only when the Bloch wave number is sufficiently small, so that g/( « 1, 
where h is the chain period, can we neglect the mathematical distinction between the Bloch and the sinusoidal waves 
and treat an LPC as an essentially continuous system. How strong the above inequality should be is not evident a 
, priori; we will comment on this question in the concluding part of the paper. 

Dispersion curves have been previously computed for polarization waves propagating in periodic chains of Drudean 
spherical nanoparticlesi* 2 -'^. It was found that there exist excitations with frequencies w and Bloch wave numbers 
q such that uo < qc^, Ch being the speed of light in the surrounding (host) medium. The phase velocity of such 
excitations, v p = u)/q, is less in magnitude than c/j. SPPs with \v p \ < Ch are considered to be outside of the "light 
• • ■ cone" and can propagate along the chain without radiative losses, similarly to plane waves in homogeneous dielectrics. 
We will reserve the term "SPP" specifically for this type of excitations. In a chain of perfectly conducting (lossless) 
particles, an SPP can propagate infinitely without decay. Of course, absorptive (Ohmic) losses in realistic metals 
' always result in exponential spatial decay of SPPs. 

The property of the dispersion curve cu(q) which is specific to the case of LPCs made of spherical particles is that it 
is very fla'bi'^, with only a weak dependence of the frequency on the Bloch wave number. The dispersion curves are 
particularly flat for SPPs polarized transversely to the chain. This results in very small group velocities, v g « cj,, A 
factor Vg/ch ~ 10~ 2 is typical. One practically important consequence of the dispersion curve flatness is a relatively 
narrow SPP bandwidth. That is, an SPP can be excited in a chain by a spatially-localized external source only in a 
narrow band of frequencies. This can be expected to significantly limit the potential application of LPCs as optical 
waveguides. 

Maier et al. have pointed out that the use of spheroidal rather than of spherical nanoparticles can result in an 
increased SPP bandwidth and, correspondingly, in a longer propagation distance 12 . In the above work, the bandwidth 
was defined as twice the spectral shift of the nearly homogeneous SPP (characterized by q — 0) with respect to the 
plasmon peak of an isolated nanoparticle and it was assumed that the propagation distance is proportional to the 
inverse of the latter. The results were confirmed by FDTD simulations in an LPC of seven prolate nanospheroid 
whose longer axis was perpendicular to the chain. 
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In this paper, we further investigate the effects of nonsphericity of the LPC constituents. We compute the dispersion 
relations in such LPCs in the dipole approximation. We find that the dispersion curves in LPCs are dramatically 
altered by replacing spherical particles with prolate or oblate spheroids. In particular, the SPP bandwidth can be 
significantly increased, in agreement with the results of Maier et at. Here, however, we define the bandwidth as the 
range of frequencies in which efficient SPP transport along the chain is possible. This includes modes with various 
values of q, including those with q ~ n/h. We also find that the use of oblate spheroids (nanodisks) whose shorter 
semiaxis is parallel to the chain is even more beneficial as it allows to achieve the desired effect at relatively modest 
values of the aspect ratio. The increased bandwidth is expected to result in a higher maximum bit rate and longer 
propagation distances for signals transported by an LPC waveguide. 

We further report that at some critical values of the spheroid aspect ratio, gaps appear in the first Brillouin zone of 
the lattice. Propagating SPPs do not exist when the Bloch wave number q is inside one of such gaps. When q is near 
the gap edge, a number of interesting phenomena take place. First, the propagation distance (the decay length) is 
dramatically increased, as compared to the same quantity when q is far from the edge. Second, the dispersion curves 
acquire very large positive or negative slopes. In this case, relatively fast (ch — \v g \ *C ch) and even superluminal 
(\v g \ > Ch) wave packet propagation can be obtained. Note that superluminal group velocity does not contradict 
special relativity^. Superluminal wave packets exist in nature and were observed experimentall y 14 ' 15 . 

Theory and numerical simulations presented below are based on the dipole approximation. The use of this ap- 
proximation dictates that the inter-particle spacings are larger than a certain threshold at which excitation of higher 
multipole moments in nanospheroids becomes non-negligible. This has limited the range of chain parameters consid- 
ered in this paper. We, however, expect on physical grounds that using chains with smaller inter-particle separations 
may be beneficial. For instance, we expect that the SPP propagation length in such chains can be increased. However, 
in order to obtain quantitative results in that limit, a considerably more complex mathematical formalism must be 
used. The latter has been developed by Park and Stroud 5 - for chains of spherical nanoparticles in the quasistatic limit. 
Unfortunately, generalization of this formalism to non-spherical particles and beyond the quasistatic approximation 
(which we deem to be essential for describing SPP propagation in long chains, as was confirmed recently in the 
experiments by Koenderink et al^) appears to be problematic. 

The paper is organized as follows. In the Section[TTl we describe and justify the basic model used to simulate Bloch 
waves and wave packets in LPCs. In Section IIII1 we compute the dispersion curves for SPPs in chains of prolate 
and oblate nanospheroids. In Section ITVl we discuss the attenuation of SPPs due to Ohmic losses in LPCs and, for 
comparison, in metallic nanowires. In Section [Vl we describe direct numerical simulations of wave packet propagation 
in LPCs. Section [VII contains a summary and a discussion of obtained results. 

II. THE BASIC MODEL 

Consider a linear periodic chain of identical metallic spheroids with semiaxes a and b (a > b). We will discuss below 
two different cases. In the first case, the chain is made of prolate spheroids whose axis of symmetry (which coincides 
with the longer axis) is perpendicular to the chain. In the second case, the chain is made of oblate spheroids whose 
axis of symmetry (which coincides with the shorter axis) is parallel to the chain. In both cases, the longer axes of the 
spheroids are perpendicular to the chain and the eccentricity, e, is given by 

e = v /i-(V«) 2 - (i) 

We will refer to the ratio b/a < 1 of the shorter and longer semiaxes of the spheroids as the aspect ratio. 

The spheroids are centered at the points x n — hn, where n is an integer and h is the chain period. The surface- 
to-surface separation of two neighboring spheroids is a = h — 2b; we require that h > 2b to avoid geometrical 
intersection of particles. Note that a stronger condition on the interparticle separation will be imposed below. The 
smaller semiaxis b is assumed to be on the order of lOnm while a can be up to a few times larger. We will see that 
SPPs propagating in such chains have frequencies to such that the corresponding wavelength in the host medium, 
A = 2nch/<-<J, is considerably larger than both a and b. Accordingly, we adopt the dipole approximation. 

This choice, however, requires an additional justification. While the dipole approximation accuracy for electro- 
magnetically interacting spheroids has not been studied directly, many results are available for spheres. The most 
basic and frequently considered example is that of two electromagnetically-interacting spheres in close proximity of 
each othe r 17 ' 18 . In this case, the spatially inhomogeneous fields scattered by the spheres result in the excitation of 
vector spherical harmonics of all orders, even if the spheres are small compared to the wavelength. More specifically, 
the electric field inside each sphere can be expanded into the vector spherical harmonics with nonzero coefficients 
appearing in arbitrarily high orders. In the dipole approximation, only the first-order terms (I = 1, m = 0, ±1) are 
retained in this expansion. The accuracy of this approximation was found to be dramatically affected by polarization. 
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If the electric field polarization is parallel to the axis connecting the spheres, the dipole approximation starts to 
deviate from the exact solution when a as Q.5R, R being the sphere radius, for both dielectric^ and conducting^ 
spheres. However, if the polarization is perpendicular to the axis, the dipole approximation yields results (i.e., the 
total dipole moment of the spheres^) with a relative error of only 2% even in the case a = 0. A careful studyi^ of 
the transversely-polarized electromagnetic modes of finite-length linear chains of interacting spheres has revealed that 
the effect of multipole interaction is to slightly shift and broaden the dipole resonance - an effect hardly observable 
in most materials due to the spectral line broadening associated with Ohmic losses. 

Below, we work in the regime when a = 2b (h = 46). In the case of spheres (6 = R), one could expect the dipole 
approximation to be very accurate for such relative separations regardless of polarization. We can, however, apply a 
more stringent test and compare a to the larger semiaxis of the spheroid, a. Within the dipole approximation, the 
physical effect described below is manifest for the following aspect ratios. For transversely polarized SPP, the gap 
in the first Brillouin zone of the chain lattice appears when b/a < 0.25 in the case of prolate spheroids and when 
b/a < 0.35 in the case of oblate spheroids. For longitudinal SPP polarization, the gap appears when b/a < 0.1 in the 
case of prolate spheroids and when b/a < 0.25 in the case of oblate spheroids. 

Consider first the transverse polarization. The aspect ratio b/a = 0.25 corresponds to a/a = 0.5. In the case of 
two spheres of radius R separated by the surface-to-surface distance a = 0.5i?, the multipole effects are negligible. In 
fact, even a smaller aspect ratio of b/a — 0.15, which is used in the numerical simulations of transversely-polarized 
wave packets in Section [Vj corresponds to the relative separation a = 0.3R. At this separation, the multipole effects 
can still be safely ignored. In the case of the longitudinal polarization, the aspect ratio which is required to observe 
the effect in chains of prolate spheroids is 0.1 which corresponds to a = 0.2i?. The multipole effects in this situation 
are expected to be significant but not dramatic. However, in oblate spheroid chains, the required aspect ratio is 0.25 
which corresponds to a = 0.5i?. At this relative separation, the effects of higher multipoles are noticeable but small, 
even for the longitudinal polarization. 

Finally, a simple physical explanation for the dramatic polarization dependence of the dipole approximation accuracy 
is available. When the polarization is parallel to a chain of spheres, the sphere surfaces which are adjacent to the 
"junctions" are similar to usual capacitors and acquire large and sign-opposite surface charge densities which, in turn, 
results in highly non-uniform, strongly enhanced local fields. This causes excitation of very high multipole moments. 
However, for the case of transverse polarization, the surface charge densities near the junctions are proportional to 
the geometrical factor cos 9 (9 being the angle between the polarization vector and the radius vector of a point on the 
sphere surface) and are small. Obviously, this consideration holds for spheroids as well. 

We thus conclude that the use of the dipole approximation is well justified for the purpose of this paper. In the case 
of transverse SPP polarization, the approximation accuracy is exceedingly good. For the longitudinal polarization, 
the accuracy can be questioned in prolate spheroid chains but is quite reasonable when oblate spheroids are used. We 
will report numerical computation of the dispersion curves for both prolate and oblate spheroids, in transverse and 
longitudinal polarization, and for various aspect ratios of the spheroids ( Section UTTT Figs. I2I4I3I5[) . However, direct 
simulation of wave packet propagation (Sections|Vl Figs. 1819ft is reported only for the choice of parameters such that 
the accuracy of the dipole approximation is not in doubt. 

In the dipole approximation, each nanoparticle is characterized by a (possibly, tensor) dipole polarizability a(w) 
and radiates as a point dipole. The Cartesian components of the nanoparticle dipole moments d n are coupled to each 
other and to the external electric field by the coupled-dipole equation£ ' 20 ' 21 i 22 , which we write here in the frequency 
domain as 



Here ££ xt is the external field amplitude at the n-th site, k = oj/ch is the wave vector in the host material at the 
frequency to and Gk{x, x') is the appropriate element of the frequency-domain free-space Green's tensor. In this paper, 
we consider both the transverse and the longitudinal polarizations of the SPP. In the absence of magnetic polarizability 
of the nanoparticles (which is assumed), the SPPs with the three orthogonal polarizations are not electromagnetically 
coupled to each other. Therefore, each polarization can be considered separately and the quantities appearing in 
Eq. $2$ should be understood as follows: d n and £^ xt are projections of the dipole moments and of the external 
electric field on the selected polarization axis, u(ui) is the appropriate scalar element of the polarizability tensor and 
Gk(x,x') is defined by 




(2) 
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- exp(ik\x — x'\) , transverse polarization , 
- , 3 ) exp(ik\x — x'\) , longitudinal polarization . 
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III. THE DISPERSION RELATIONS 

An SPP mode is an excitation that propagates along the chain without an external source. Thus, to find the 
dispersion relation, we seek a solution to ^ with zero free term, -E® xt = 0, in the form d n cx exp(iqx n ), where q is in 
the first Brillouin zone of the lattice, q € [— Tr/h, ir/h]. Substitution of this ansatz into ([2]) yields the equation 



oT 1 {uj) = h~ 3 S(hk,hq) , 
where S{hk, hq) is the dimensionless dipole sum (the dipole self-energy) of the chain defined by 
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exp(m£) cos(nry) , transverse polarization , 
exp(in^) cos(nrj) , longitudinal polarization . 



(5) 



Note that the dipole sum is a function of two dimensionless parameters £ — kh and r\ = qh but does not depend on 
the particle shape and material properties. 

The dispersion relation, i.e., the mathematical dependence of the SPP frequency u> — kch on its Bloch wave number 
q, can be obtained by finding all pairs of variables (u>,q) that satisfy Eq. (j4|). This can give rise to one or more 
branches of the complex function Lu(q). Generally, purely real pairs (w,g) that solve ((J) do not exist. One can, 
however, consider purely real values of q and seek complex frequencies, as was done by Koenderink and Polman 3 . The 
imaginary part of lu is then interpreted as the SPP decay rate. Alternative approaches include numerical computation 
of discrete modes in a finite chain 1 and plotting Im[a _1 — h^^S]^ 1 as a function of two variables k and q and visually 
identifying the points at which this function appears to have a maximum or a saddle point^. 

In this paper, we are interested in propagation of wave packets which are excited as superpositions of oscillations 
with purely real frequencies. Since SPPs propagate without radiative losses, their Bloch wave numbers q arc purely 
real if the chain is made of a non-absorbing material such as an ideal conductor. However, when Ohmic losses in 
realistic metal are accounted for, q acquires an imaginary part. In what follows, we use two different approaches to 
computing the dispersion curve. In Section IIII A[ we consider chains made of ideal (lossless) metal and seek purely 
real solutions w(q), as was suggested by Simovski 2 . Numerically, this is accomplished by finding pairs of real variables 
{co,q) that satisfy the dispersion equation ^ by the method of bisection. Such purely real solutions exist if the 
permittivity of metal is taken to be real. Next, in Section [ill A[ we consider realistic metals. Here we seek pairs (w, q) 
that satisfy the dispersion equation such that oj is purely real but q is complex. Numerically, such pairs are obtained 
by utilizing the root-finding algorithm implemented in Wolfram's Mathematica. For the specific case of silver, we 
find that both approaches yield the results which are very close quantitatively when the dependence of u> on Re(g) is 
considered; the first approach, however, provides no information on SPP attenuation which is governed by Img. In 
both cases, we seek solutions only in the region Reg > k = Lo/ch\ as was discussed in the Introduction, excitations 
with Reg < k experience radiative decay in addition to Ohmic losses and are not considered in this paper. 

A. Dispersion Relations for Ideal Metal 

In this Section we assume that q and k are real and view the dipole sum S(kh, qh) as a function of two purely real 
variables. As the first step, we write the inverse polarizability in (T?|) in the form—: 
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where q;ll (w) is the Lorentz-Lorentz quasistatic polarizability of nanoparticles and 2ifc 3 /3 = i(2/3)(uj/ch) 3 is the first 
non-vanishing radiative correction to the inverse polarizability. The latter is given by 



-1/ \ 47r ( , e h \ 

where v is the appropriate depolarization factor and v is the spheroid volume^. 
In the case of prolate spheroids, the volume is given by 



4tt 



-ab A 



(8) 



and the three depolarization factors are v\ for polarization along the spheroid axis of symmetry and V2 = vz = (1— fi)/2 
for two linearly independent transverse polarizations, where 



1-e 2 



1 . 1 + e 
— In 

2e 1-e 



(9) 



Here e m and are the permittivities of the (metallic) spheroids and of the host medium, respectively, and e is the 
spheroid eccentricity given by formula ([Tj). 
For oblate spheroids, the volume is 



and the depolarization factors are v\ — V2 for two linearly independent polarizations which are orthogonal to the 
spheroid axis of symmetry and 1/3 = 1 — 1v\ for the polarization along the axis of symmetry, where 



5(e) 
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(11) 



While both permittivities e m and eh have, in general, some frequency dependence, here we neglect the dispersion in 
the host and assume that eh = const > 0. Then, for nanoparticles made of a lossless material, Im(a^) — 0. At the 
same time, if q > k, the imaginary part of the dipole sum is&2£ lm[S(kh, qh)] — —2(kh) 3 /3 and, in the region of (fc, q) 
which is of interest to us, Im [a -1 — h~ 3 S~\ — 0. Therefore, the imaginary part of Eq. ^ is satisfied identically and 
only its real part needs to be considered. We then utilize Eqs. ©,® and arrive at the following dispersion equation: 



v + Re ( 6h ) = e h -^Re[S(hk, hq)} . (12) 

In the case of a lossless metal and a transparent host medium, the real-part symbol in the left-hand side of (fT2"]) can 
be omitted. 

Eq. (fT2|) allows one to analyze the relation between the spheroid aspect ratio and the dispersive properties of the 
LPCs. Consider first the SPP polarization which is transverse to the chain. For this polarization, the dependence of 
the real part of the dipole sum S(kh, qh) on its arguments is illustrated in Fig. [lja) and the relevant depolarization 
factor is v\ (for both prolate and oblate spheroids). We now notice the following. In the spectral region where the 
metal experiences anomalous dispersion, the second term in the left-hand side of (fT2]) is negative, assuming that the 
host medium is a transparent dielectric. When the aspect ratio b/a is decreased, the depolarization factors v\ given 
by Eqs. (0 or (fTT|) . for prolate and oblate spheroids, respectively, both approach zero. As a result, the whole left-hand 
side in (|12[) becomes negative. On the other hand, there are regions in the (fc, q) space in which the right-hand side 
of (fT2")l is strictly positive. Thus, it can be seen from Fig. [Ija) that Re[S(kh, qh)] > if qh/ir > 0.5 and q > k. As 
a result, for sufficiently small ratio b/a, equation (|12p ceases to have real- valued solutions if q > q^r, where q^r is the 
(aspect ratio-dependent) critical value of q specific to the transverse polarization. The interval of Bloch wave numbers 
q > q^r corresponds to a gap in the first Brillouin zone of the chain lattice in which SPPs do not exist. It is important 
to emphasize that the critical constant q^r < ir/h exists only for sufficiently small ratio b/a. 
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FIG. 1: (color online) Contour plot of the real part of the dipole sum, Ke[S(kh, qh)] for the transverse (top) and the longitudinal 
(bottom) polarizations. Due to the symmetry S(£,ri) = S(£, — rj), only the positive half of the first Brillouin zone is shown. 
Note that, for the transverse polarization, Ke[S(kh, qh)] diverges logarithmically (approaches positive infinity) on the light line 
q — k. Near this line, the function changes so fast that it is not feasible to depict it quantitatively using the contour plot; the 
region of divergence is schematically shown as a diagonal white line in the top panel. There is no divergence for the longitudinal 
polarization. 

Similar considerations can be applied to the longitudinal SPP polarization, except that the relevant depolarization 
coefficient is, in this case, z/3. The dependence of Ke[S(kh, qh)] on its arguments is illustrated in Fig. [ljb). It can 
be seen that Ke[S(kh, qh)] is strictly positive for qh/n < 0.4 and q > k. Correspondingly, Eq. (fT2")) ceases to have 

real- valued solutions for q < q| which defines a gap in the first Brillouin zone of the lattice. Here ql is the critical 
Bloch wave number for the parallel polarization; as in the case of transverse polarization, gj| exists only for sufficiently 
small aspect ratios and is aspect ratio-dependent. Note that the depolarization factor V3 approaches a finite value 
rather than zero when the aspect ratio is decreased. This limit is 1/2 for prolate and 1 for oblate spheroids. Due to 
this reason, observation of the gap for longitudinally-polarized SPPs requires a smaller aspect ratio. This point will 
be illustrated below by numerical examples. 
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We now compute the dispersion curves numerically. To solve Eq. (|12p . a specific expression for the metal permittivity 
e m is needed. We use here the Drude formula 



£ m = e o -, — — -r , (13) 

UJ[bJ + IJ) 

where eo is the contribution due to interzone transitions^, u p is the plasma frequency, and 7 is the relaxation constant. 
In this subsection, we set 7 = to describe a lossless metal (realistic values of 7 will be used in Sections IIIIBI and 
IVl below) . We then solve Eq. (fT2|) by the method of bisection to obtain the dispersion curves and the SPP group 
velocity. 

The results are shown in Figs. [21 through [5j The top panels in these four figures show a number of dispersion 
curves computed for different aspect ratios b/a, as labeled, and for different SPP polarizations. The group velocity 
of the SPPs, v g — duj/dq, is plotted as a function of q in the bottom panels. The values of eh and eo are taken to be 
2.5 (as in the case of a glassy medium) and 5.0 (the experimental value for silver—), respectively. The computation 
does not depend on the absolute value of the plasma frequency uj p but on the dimensionless parameter X p /h, where 
X p = 2izCh/u) p is the wavelength (in the host medium) at the plasma frequency. This ratio was chosen to be X p /h = 3.4. 
Finally, we have set the ratio h/b = 4 so that the minimum surface-to-surface separation of two neighboring spheroids, 
a, was equal to 2b. As an illustration, for the specific case of silver, we have the following parameters: the vacuum 
plasma wavelength is Ap™ "* « 136nm and the corresponding value in the host medium is A p = ' / \fth ~ 86nm; 
correspondingly, h ~ 25nm, b « 6nm and a varies from 6nm to 60nm. The latter value was obtained for the smallest 
aspect ratio used, b/a = 0.1. It can be seen that, for all points on the dispersion curves shown in Figs. SPP 
frequencies are well below the plasma frequency to p and kh/n <C 1. The last inequality is especially strong for 
the central frequencies (indicated by horizontal arrows in Figs. [2] and [5]) which were used to simulate wave packet 
propagation in Section [V] below. This reconfirms the dipole approximation validity. 

We now discuss the computed dispersion curves in more detail, starting with the case of transverse SPP polarization 
(Figs. HE]). First, we note that in infinite, strictly periodic chains, the dispersion curves of transversely-polarized SPPs 
start at the point k = q = and then follow the light line k — q for some range of q. This small-g part of the dispersion 
curve is related to the logarithmic divergence of the dipole sum on the light line^ and is extremely difficult to find 
numerically. Indeed, in order to satisfy Eq. (|12p . the points on the small-q section of the dispersion curve must be 
specified with exponentially large numerical precision. This is why the small-g section of the dispersion curve has not 
been reported in a number of numerical investigations^'^. However, the SPPs with q « k -C n/h exist and were 
observed in numerical simulations^. 

In this paper, we do not consider the small-q part of the dispersion curve but focus on the SPPs for which the ratio 
k/q is considerably less than unity. Such modes exist for qh/ir > 0.2 for all four values of b/a shown in Figs. 12131 The 
main point of this paper is that the dispersion curve shape is strongly influenced by the ratio b/a. When b/a = 1.0, 
the corresponding dispersion curve is almost flat (apart from the linear small-g section of the curve). At b/a = 0.5, 
the curve begins to bend down noticeably at larger values of q. Finally, when b/a < 0.25 (in the case of prolate 
spheroids) or when b/a < 0.35 (in the case of oblate spheroids), the curves cross the k — axis at the point q = q^ 
and no solutions exist for q > qj: . Near the critical point, the dispersion curves acquires a very large negative slope. 
Note that the corresponding slope is positive in the q < part of the first Brillouin zone (which is not shown in the 
figures) . 

From comparison of Figs. 12131 a conclusion can be made that the dispersion curves are more sensitive to the aspect 
ratio in the case of oblate spheroids. In particular, the gap in the first Brillouin zone of the chain lattice appears for 
more moderate values of b/a if the chain is made of oblate spheroids. 

Dispersion curves for the longitudinal SPP polarization are shown in Figs. 14151 It can be seen that, for sufficiently 
small aspect ratios, there exists a critical value qt such that Eq. (| 1 2[) has no real- valued solutions for q < q\ . Thus, 



SPP propagation in the chain is only possible for q larger than the critical value q\ (if the latter exists). The group 

velocity shown in Figs. IHb)J5]Jb) acquires large positive values in the vicinity of q\. 

Similarly to the case of transverse SPP polarization, the dispersion curves are more sensitive to the aspect ratio if 
the chains are made of oblate spheroids. Thus, in the case of prolate spheroids, the gap in the first Brillouin zone of 
the lattice appears only when b/a < 0.1. However, if the chain is composed of oblate spheroids, the gap appears when 
6/a<0.25. 
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FIG. 2: (color online) Dispersion curves (a) and group velocity (b) for transversely polarized SPPs in chains built from prolate 
spheroids whose axis of symmetry is perpendicular to the chain, for different spheroid aspect ratios b/a and for fixed ratios 
h/b = 4 and X p /h = 3.4. Since the dispersion curves are symmetric with respect to the fc-axis, only the positive half of the 
first Brillouin zone (q > 0) is shown. Only data points that were found numerically are plotted. Theoretically, however, all 
dispersion curves in infinite chains start at the point k = q = and follow the light line (labeled as k — q in the top panel) for 
some range of q's. Horizontal arrows indicate central frequencies of the Gaussian wave packets whose simulated propagation is 
illustrated in Fig. [8] below. 



B. Dispersion Relations for Realistic Metal 

The dispersion curves shown in Figs. [2] through [5] contain no information on either the rate or the direction of SPP 
spatial decay. However, it can be seen in these figures that an SPP is characterized by phase and group velocities 
and that these can have different signs when projected onto the x-axis. On physical grounds, we expect that wave 
packets should decay in the direction of propagation. This imposes certain restrictions on the signs of the real and 
the imaginary parts of q. Thus, if v g v p < 0, we expect that Ke(q)lm(q) < while if v g v p > 0, we expect that 
Re(q)lm(q) > 0. 

The above statement can be illustrated by considering the following thought experiment. Assume that a short 
Gaussian optical pulse is injected into the central part of a long chain by a spatially localized source such as a near- 
field microscope tip operating in the illumination mode. The pulse will propagate in the form of two wave packets 
in both directions along the chain. If v g v p < 0, both wave packets would be composed of Bloch waves whose phase 
velocities point towards the source and group velocities point away from the source. Since the sign of the phase 
velocity is the same as that of Re(g), and since both wave packets should decay in the direction of propagation, we 
expect in this case that Re(g)Im(g) < 0. Similar consideration can be applied to the case v g v p > 0. 

In the numerical simulations of this subsection, we have verified that the above analysis is, indeed, correct. To this 
end, we have incorporated the Ohmic losses into the model of metal permittivity. Specifically, we have set the Drude 
relaxation constant in Eq. (fT3|) to be 7 = 0.002^ which is the experimental value for silver. Note that Eq. (| 1 2[) can 
still be satisfied in this case with a purely real pair of (to, q) (the real-part symbol in the left-hand side of the equation 
must be retained if e m is complex- valued) . However, the imaginary part of the more general Eq. Q is no longer 
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FIG. 3: (color online) Same as in Fig. [2] but for a chain made of oblate spheroids whose axis of symmetry is parallel to the 
chain and for a different set of aspect ratios. SPP polarization is orthogonal to the chain. 

satisfied identically. We, therefore, must seek such pairs of variables (w, q), where q is now complex, that satisfy both 
the real and the imaginary parts of Eq. Q . This can not be achieved by the use of the simple bisection algorithm 
that was employed in the previous subsection. Instead, we employ here the Wolfram's Mathematica root finder to 
obtain complex roots of Eq. (j4]), q, for each real value of u>. 

The results of this simulation are plotted parametrically in the complex g-plane in Fig. [5] for a chain of prolate 
spheroids of the aspect ratio b/a = 0.15 and in Fig. [7]for a chain of oblate spheroids of the aspect ratio b/a = 0.25. 
Other parameters are the same as in Fig. [3J As expected, the real and imaginary parts of the Bloch wave number 
have opposite signs for the transversely-polarized SPPs so that the wave packets decay in the direction of propagation 
specified by v g , even though the phase velocity is pointing into the opposite direction. For longitudinally polarized 
SPPs the product v p v g is always positive and, correspondingly, Im(q)Re(g) > 0. 

An important observation is that when Keq approaches one of its critical values, Imq tends to zero. Correspondingly, 
the decay length of the SPPs is dramatically increased. 

We note that the data shown in Fig. [6] represent essentially the same dispersion curves as the ones plotted in 
Fig. [2ja) and[4ja) for the case b/a = 0.15. Similarly, data in Fig. [7] correspond to the dispersion curves plotted in 
Fig. GHa) and[^a) for b/a — 0.25. The discrepancy between the curves <x>[Re(g)] computed by the two methods is 
negligibly small. This further justifies the utility of the simple numerical method of Section IlII Al for computing the 
dispersion curves in LPCs. 

IV. ATTENUATION OF SPPS AND THE DECAY LENGTH 

The nonzero relaxation constant in Eq. (| 13[) leads to exponential decay of SPPs even in the absence of radiative 
losses. The decay lengths in LPCs is given by £ LPC = l/lm(q). If q is not very close to the edge of one of the gaps that 
were discussed above, £ hPC can be computed in the quasi-particle pole approximation^ which results in the following 
expression: 
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FIG. 4: Same as in Fig. [2] but for longitudinal SPP polarization and a different set of aspect ratios. Note that, unlike in the 
case of transverse polarization, the dispersion curves do not have linear small-q segments. As stated in the text, only data 
points that satisfy q > k are plotted, since there are no real-valued solutions in the region q < k. 



dReS(£, rj) 



drj 



(14) 



where <!?(£, 77) is the dipole sum given by Eq. ([5]) and qo(£) is a purely real solution to (JT2J) for a given value of the 
dimensionless parameter £ — toh/ch- Under the assumptions that e/j,Im(e m ) <C Re(e m ), we have 



- Im KiM ~ ■ (15) 

As we have seen above, the real-valued solutions to (I12|) satisfy, albeit approximately, the more general dispersion 
equation (H}. Therefore, we can evaluate the derivative in Eq. (fT4"l) at the point (£,77) with the understanding that £ 
and r\ are purely real variables that satisfy the approximate dispersion equation \Y1\ . We further neglect the terms 
that are of the order of ~ l/£ and ~ l/£ 2 in the square brackets of the expression ([5]) for the dipole sum and arrive 
at the following estimate: 



, Av 
2ttIi 3 luj 



cos(ujhm/ch) sin(qhm) 



(16) 



where A = 1 for the transverse polarization and and A = 2 for the longitudinal polarization. 

A direct calculation for the transversely-polarized SPP in a prolate spheroid chain, b/a = 0.15 and other parameters 



same as in Fig. [5] yields the decay lengths £ LPC ~ 7/jm for u 



O.ltjp and 



15/jm for uj — 0.05w p . Interestingly, 



the polarization dependence of decay length is contained solely in the factor A. However, the two polarizations have 
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FIG. 5: Same as in Fig. [3] but for longitudinal SPP polarization and a different set of aspect ratios. Note that, unlike in the 
case of transverse polarization, the dispersion curves do not have linear small-g segments. As stated in the text, only data 
points that satisfy q > k are plotted, since there are no real-valued solutions in the region q < k. Horizontal arrows indicate 
central frequencies of the Gaussian wave packets whose simulated propagation is illustrated in Fig. [9] below. 



markedly different dispersion curves and therefore, same pairs of variables (£, rf) may not be accessible for the two 
different polarizations. 

It is instructive to compare the SPP decay length in nanoparticle chains and metallic nanowires. The dispersion 
equation in a metallic cylindrical waveguide i o 27 ' 28 : 



€ m h(K m R) + e h Ki(K h R) _ Q 

(K m R) K h K (K h R) 

where R is the cylinder radius, h{x) and Ki(x) are the modified Bessel functions, K mt h = \/<Z 2 ^m,h^> 2 /c| and the 
indices "m" and "/i" label the quantities for the metal and for the surrounding dielectric host, respectively. If the 
wire is sufficiently thin, we can expand the Bessel function to the first non- vanishing order in n m R and n^R to obtain 
the simplified dispersion equation 



1 ( 

6m {n h Rf \n{n h R/2) + C ' 1 ' 

where C is the Euler constant. This equation can be solved approximately (with logarithmic precision) as 
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FIG. 6: Dispersion curves plotted parametrically in the complex g-plane for transversely (a) and longitudinally (b) polarized 
SPPs in an LPC of prolate spheroids. Parameters: 'y/ujp — 0.002, b/a — 0.15 and other parameters same as in Fig. [2] Dots 
label the values of the dimensionless parameter kh/n = uih/irch- In panel (a), only the points that correspond to the the region 
VpVg < are shown, which corresponds, approximately, to \Re(q)\h/iv > 0.17. 



We then use the Drude formula for e m , take into account the fact that the propagation constant Reg in a metal 
nanowire is much larger than ^JThUj / c (which is the wavenumber in the surrounding medium) and obtain the following 
estimate for the decay length: 



(20) 



One additional condition that has been used in deriving the above estimate is u> 3> 7. For a silver nanowire of radius 
R = 25nm, the estimate yields £ w i rc ~ 2.8/im. 

The following conclusions can be made. While the SPP decay length in nanowires depends only on the metal 
and host permittivities and the waveguide radius, the same quantity in the LPCs can be controlled by changing 
the inter-particle separation and the particle dimensions. Decreasing the inter-particle spacings results in stronger 
electromagnetic coupling and larger propagation distances. However, for sufficiently small values of h, the dipole 
approximation breaks down and the estimate (|TB)) becomes invalid. The dimensions of a nanoparticle provide another 
set of degrees of freedom in controlling the decay length in LPCs. Note that the expression (fTB)) contains an overall 
factor v/2irh 3 which can be interpreted as the volume fraction of metal (in a unit cell of the chain lattice). For an 
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FIG. 7: Same as in Fig. [6]for an LPC of oblate spheroids whose aspect ratio is b/a = 0.25 

LPC of prolate spheroids, this factor is equal to 2ab 2 /3h 3 . However, for an LPC made of oblate spheroids, the factor 
is 2a 2 b/3h 3 . It obtains that the propagation distance in oblate spheroid LPCs is effectively increased by the factor of 
a/b > 1 compared to prolate spheroid LPCs. Thus constructing an LPC from thin nano-disks whose axis of symmetry 
is parallel to the chain axis may be beneficial. 

V. MODELING OF WAVE PACKET PROPAGATION 

We now demonstrate that superluminal wave packets can indeed propagate in chains with sufficiently small aspect 
ratios b/a. To this end, we consider a finite chain of N nanoparticles excited by a pulse with Gaussian temporal 
profile incident on the first particle in the chain. In the time domain, the pulse is described by the formula 



E^(t) = S nl £ exp [-iuot - (t/At) 2 ] , 
where £ is an arbitrary amplitude and At is the pulse duration. This function has the Fourier transform 



(21) 



E%?(w) = <5„i^At£exp 



At 



(22) 
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FIG. 8: (color online) Envelopes of transversely polarized wave packets in a chain of N = 5000 prolate nanospheroids with the 
aspect ratio b/a = 0.15 at different moments of time t. Spheroids are oriented so that their axes of symmetry are perpendicular 
to the chain axis. Time is measured in the units of r = h/cu- Column (a): ujq = 0.1ui p , v g ~ 0.57c/j. Column (b): u>o ~ 0.05oj p , 
v g ~ 1.16cfe. Arrows indicate that the wave packet propagates from right to left after being reflected from the far end of the 
chain. 

The numerical procedure is as follows. The above expression for £ , ° xt (w) is used as the free term in the right-hand 
side of Eq. ([2]). The equation is solved numerically by direct matrix inversion for multiple values of uj sampled in 
a sufficiently large interval and with sufficiently small step to ensure convergence. This yields a family of numerical 
solutions d n (oj). The real time quantities d n (t) are then obtained by the inverse Fourier transform according to 



d n (t) = / d n (uj) exp(-iut) 



2^ 



(23) 



Numerically, this integral is evaluated by the trapezoidal rule. 

We have simulated transversely polarized wave packets in LPCs of prolate spheroids with the aspect ratio b/a — 
0.15. Longitudinally polarized SPPs were simulated in chains of oblate spheroids with the aspect ratio b/a = 0.25. 
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FIG. 9: Envelopes of transversely polarized wave packets in a chain of N = 5000 oblate nanospheroids with the aspect ratio 
b/a — 0.25 at different moments of time t. Spheroids are oriented so that their axes of symmetry coincide with the chain axis. 
Time is measured in the units of t = h/c^. Note that the largest time shown in this figure is twice smaller than the respective 
quantity in Fig. EI Column (a): ljq = 0.15oj p , v g m 0.77ch. Column (b): u>o = 0.05oj p , v g m 2.57ch- Arrows indicate that the 
wave packet propagates from right to left after being reflected from the far end of the chain. 



Simulations were performed in a chain consisting of N = 5 ■ 10 3 spheroids; the overall length of the chain was (given 
h = 25nm) L — 125/im. All parameters were the same as those used for calculating dispersion curves shown in Figs.[2l 
with the only exception that Ohmic losses in the metal were taken into account by means of using the nonzero Drude 
relaxation constant 7 = 0.002cj p . Four sets of simulations have been performed, the first two for the transverse and 
the other two for the longitudinal polarization. 

In the case of the transverse polarization, two different central frequencies of the pulse have been used. The first 
pulse had the central frequency luq = O.lujp (correspondingly, koh/n — 0.06 where fco = <-^o/ch) and the pulse spectral 
width was Au> = ujq/5. Thus the excitation was relatively broad-band but very narrow in the time domain: given the 
experimental value of u> p for silver, the above spectral width corresponds to At sa 7.2fsec. The second pulse had the 
central frequency twice smaller than the first, loq = 0.05w p , with the same relative spectral width Aw = u>o/5. In the 
time domain, this corresponds to At » 14.2fsec. The central frequencies of the two pulses are shown in Fig. [5](a) by 
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the horizontal arrows. Amplitudes of the obtained wave packets are illustrated in Fig. [8] at different moments of time 
measured in the units of t = h/c^. The maximum time shown on the plots is t — 6 10 3 r 800fsec. 

It can be seen that the wave packet with loq = 0.1w p propagates away from the source with the subluminal group 
velocity v g ~ Q.blch- However, the wave packet with the smaller central frequency propagates at the speed v g ~ 1.16c/j. 
These group velocities are in quantitative agreement with the data shown in Fig. [2^b) . Note that the group velocities 
can be evaluated as the slopes of the dispersion curve drawn for b/a = 0.15 in Fig. [2a) at the central frequencies 
indicated by the horizontal arrows. As expected, the superluminal wave packet is spreading faster than the subluminal 
wave packet. This is so because of the larger value of the second derivative d 2 uj/d 2 q at the smaller central frequency. 
Yet, near the end of the chain, the duration of the superluminal pulse is still only « lpsec. 

The two sets of simulations for the longitudinal polarization are shown in Fig. [5J Here we have used a chain of 
oblate spheroids with the aspect ratio b/a = 0.25. The central frequencies of the two pulses were uj = 0.25cl> p and 
luq = O.lujp. The pulses' relative spectral widths were the same as in the case of the transverse polarization, namely, 
the pulse spectral width was Alu — luq/5. By tracing the maximum of each wave packet, we deduce v g ~ 0.88cft, for 
ujo = 0.25w p and v g = I.YJch for ujq = O.lujp. This is in full agreement with the group velocities shown in Fig.[5jb). 

Finally, note that the decay lengths in each simulation can not be easily deduced from the time evolution of the 
maxima of the wave packets. This is because the propagation is accompanied by both decay and spreading. The 
latter takes place even in the absence of Ohmic losses. 

VI. CONCLUDING REMARKS 

In this paper, we have employed the coupled dipole approximation to compute the dispersion curves and to 
model propagation of wave packets of surface plasmon polaritons (SPPs) in linear periodic chains (LPCs) of metallic 
nanospheroids. The main novel element of this study, as compared to the previous work on the subject 1 -'^' 3 -' 4 -' 5 -'^'^, is 
the account of particle nonsphericity. 

We have shown that the group velocity, decay length and the bandwidth of surface plasmon polaritons (SPPs) 
propagating in linear periodic chains (LPCs) of metallic nanoparticles can be effectively tuned. The tunability is 
achieved by means of varying the nanoparticles aspect ratio. The decay length can be dramatically increased for 
Bloch wave numbers q near the edges of gaps that appear in the first Brillouin zone of the lattice for sufficiently small 
aspect ratios. At the same time, the SPP group velocity is increased up to superluminal values. By replacing the 
host medium with vacuum, it is also possible to excite a wave packet whose group velocity is larger than the speed of 
light in vacuum. Such wave packets exist in nature and were observed experimentall y 14 ! 15 . 

Comparison of Figs. [2] through [5] reveals that the parameters of two different LPCs can be tuned so that one 
LPC supports transversely polarized SPPs and the other chain supports longitudinally polarized SPPs with the same 
electromagnetic frequency. This fact can be utilized for guiding the SPPs through corners (ninety-degree turns in 
an LPC) and/or for splitting and coupling the SPPs at T-junctions. Another potentially interesting element of 
an integrated photonic circuit a two-segment straight LPC. Assume that, at a given frequency, one segment can 
support only transversely-polarized SPPs while the other segment supports only longitudinally polarized SPPs. At 
the junction of the two segments, an externally-manipulated (i.e., by magnetic field) coupling nanospheroid is placed. 
When the coupling nanospheroid makes the angle of either or 7r with respect to the chain axis, the two segments 
are decoupled and do not allow direct transmission of light pulses. However, if the coupling spheroid is rotated by the 
angle of 7r/4 with respect to the axis, the transversely-polarized SPP propagating in the first segment is coupled to 
the longitudinally-polarized SPP in the second segment and transmission along the chain becomes possible. Detailed 
investigation of these possibilities will be the subject of future work. 

In the case of transverse SPP polarization, the group and phase velocities of SPPs can be antiparallel. We, however, 
have found that the negative group velocity per se (defined here by the condition v g v p < 0) does not necessarily imply 
superluminal propagation or a negative time delay as was suggested previousl y 13 ' 29 . For example, the wave packet 
shown in Fig. [8ja) propagates slower than Ch even though it is composed of waves whose frequencies are in the 
negative dispersion region. It is important to realize that the effects theoretically described in these two references, 
as well as the corresponding experimental observations 1 ^ 5 -, involve propagation of an optical pulse from a medium 
with normal dispersion to a medium with negative dispersion and the presence of the interface is essential. In this 
paper, we are looking at a somewhat different physical situation when the optical pulse is injected into a waveguide 
by a predetermined external source which is located in the near field of the waveguide. We then observe that the 
pulse propagates away from the source with the velocity \v g \, irrespectively of the sign of the product v g v p . Thus the 
superluminal propagation is obtained when \v g \ > c^ but not necessarily when v g v p < 0. 

Antiparallel phase and group velocities that we have observed in the case of transverse SPP polarization deserve a 
separate discussion. We believe that this phenomenon can not be interpreted as "negative refraction" . The reason is 
that the LPCs considered in this paper are essentially discrete objects and can not be described by effective medium 
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parameters. The elementary excitations that propagate in LPCs are Bloch waves rather than sinusoidal waves, and 
this fact should not be disregarded. The region of negative dispersion shown in Fig. Ufa) starts at qh 0.2tt « 0.6. In 
general, the chain can be viewed as continuous only when g/(« 1. The above condition is not satisfied in the region 
of negative dispersion. It is, however, not clear a priori, how strong this inequality should be for the effective medium 
approximation to be valid. In the specific case of LPCs, one can consider the following argument. We expect the 
effective medium parameters such as the permittivity e(u>) or the refraction index n{u>) to be single- valued functions 
of their argument. However, for every point on the negative slope section of the dispersion curves shown in Fig. HJa), 
there is another point on the same curve with the same frequency but a smaller value of q. This second point is located 
on the linear, small-g section of the dispersion curve. Although this small-g section is difficult to find numerically 
(and, as a result, is often overlooked), it exists. It is therefore logical to assume that, if a chain be assigned some 
effective medium parameter for a given frequency u>, this parameter must be computed using the point on the small-q 
section of the dispersion curve. The latter exhibits positive (and linear) dispersion. In the above argument, we have 
disregarded the possibility of introducing non-local effective medium parameters which are characteristic, for example, 
of chiral media and can result in negative dispersion^ . However, the physical object that we consider in this paper 
is essentially non-chiral. 

The authors can be reached at: 
algovOseas .upenn. edu and vmarkel@mail .med.upenn. edu. 
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